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ABSTRACT 


This  volume  contains  a  detailed  description  of  the  analytical  methods 
used  in  constructing  the  Monte  Carlo  model  and  code.  The  various 
mechanisms  contributing  to  the  energy  losses  of  heavy  particles  moving 
through  the  atmosphere  are  reviewed,  and  it  is  concluded  that  for  the  cases 
of  interest  here  only  the  elastic  collisions  need  to  be  taken  into  account. 

An  outline  is  given  of  a  statistical  treatment  of  inelastic  collisions  together 
with  a  derivation  of  the  corresponding  probability  distribution.  The 
probability  distribution  for  the  beta  decays  of  fission  particles  is  also 
derived,  and  lastly,  the  logical  construction  of  the  code  for  the  calculation 
of  a  history  of  a  given  particle  is  described  in  detail. 
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INTRODUCTION 


This  report  is  a  continuation  of  a  study  carried  out  by  the  Theoretical  Physics  Group 
at  Lockheed  on  the  motion  of  the  bomb  debris.  This  volume  describes  the  mathemati¬ 
cal  background  of  the  Monte  Carlo  code  used  in  this  study  and  the  code  itself. 

We  would  like  to  acknowledge  the  contributions  of  A.  J.  Cook  to  the  formulation  and 
execution  of  the  Monte  Carlo  calculations. 
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SUMMARY 


In  any  Monte  Carlo  model  it  is  advantageous  to  limit  the  number  of  microscopic 
processes  which  have  to  be  treated  statistically.  In  die  model  of  behavior  of  fission 
particles  in  the  atmosphere,  the  most  important  processes  are  perhaps  the  collision 
of  these  particles  with  the  atmospheric  atoms.  This  is  the  subject  of  Section  1, 

Energy  Loss  of  Fissiou  Particles  Released  at  Hi^  Altitudes,  where  the  contributions 
of  the  elastic  and  inelastic  collision  to  the  slowing  down  of  a  heavy  particle  are  briefly 
examined.  It  is  found  that  approximate  calculations  of  such  contributions  are  valid 
for  velocities  up  to  about  5  x  10  cm/ sec,  and  that,  for  the  same  velocities,  a  good 
first  approximation  consists  of  neglecting  the  energy  loses  due  to  inelastic  collision. 

With  this  assumption,  it  is  possible  to  integrate  out  the  contributions  of  the  elastic 
collisions,  and  only  the  inelastic  collisions  have  to  be  treated  statistically.  The 
aiqiropriate  probability  distribution  is  derived  in  Section  2.,  Inelastic  Collisions  of 
Fission  Particles  in  the  Atmosphere.  This  distribution  depends  only  on  the  cross 
section  for  the  collision  and  the  Initial  altitude,  and  is  strongly  determined  by  the 
anal3rtical  assumption  made  on  the  behavior  of  the  density  with  altitude.  The  differ¬ 
ent  forms  of  the  distribution  in  various  atmospheres  are  briefly  examined.  The  choice 
of  the  final  form  is  governed  by  simplicity  ot  formulation  and  the  necessity  for  rapidly 
developing  a  working  code . 

In  Section  3,  Beta  Decay  of  Fission  Fragments,  a  probability  distribution  for  the 
beta  decay  of  a  typical  fission  particle  is  derived.  The  assumption  is  made  that  the 
activity  of  a  collection  of  fission  particles  is  adequately  given  by  the  Way-Wigner  law, 
and  a  prescription  for  obtaining  a  good  approximation  to  this  law  is  given. 
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In  Section  4  ,  IBM  7090Code  for  the  Monte  Carlo  Model,  the  detailed  steps  in  the  cal¬ 
culation  of  a  history  of  a  fragment  are  described.  It  must  be  pointed  out  that  all  the 
steps  shown  and  the  associated  numerical  procedures  were  chosen  primarily  to  develop 
in  a  reasonable  time  a  rapid  code.  Thus,  a  prescription  for  decay  times  differing 
slightly  from  that  of  Section  3,  was  adopted.  Section  4  also  contains  an  outline  of  the 
logical  flow  of  the  code.  Section  5  contains  a  listing  of  the  current  code  together  with 
the  output  of  a  sample  run. 
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Section  1 

ENERGY  LOSS  OF  FISSION  PARTICLES  RELEASED  AT  HIGH  ALTITUDES 


In  this  section  the  different  mechanisms  contributing  to  the  stopping  power  of  the 
atmosphere  on  neutral  and  ionized  fission  particles  are  reviewed.  The  slowing  down 
of  a  fission  particle  is  due  to  its  elastic  and  inelastic  collisions  with  the  ambient 
atoms,  and  of  the  latter  collisions  we  consider  those  giving  rise  to  ionization,  elec¬ 
tron  capture,  and  electron  loss. 

For  particle  velocities  below  about  S  x  10  cm/ sec  elastic  collisions  dominate,  the 
energy  loss  being  governed  by  the  transport  cross  section  through  the  law 


dt 


.nisi 

2 


(1.1) 


where  n(s)  is  the  number  density  of  the  ambient  medium.  The  evaluation  of 
depends  essentially  on  the  interaction  potential  between  the  colliding  systems.  For 
heavy  atoms,  the  potential  can  be  calculated  to  an  accuracy  of  10  percent  on  the  basis 
of  the  statistical  model  for  the  electrons  (Ref.  1).  Thus, 


U(r) 


21^2  « 


X{il>  r/a) 


(1.2) 


where  x(x)  is  the  screening  function  in  the  Thomas- Fermi  potential,  and 

are  the  atomic  numbers  of  the  interacting  atoms,  r  is  the  intemuclear  distance. 


a 


(97rVl28) 


1/3 


4.7  X  10  cm 


(1.3) 


4 


and 


*  -  (s 


1/2 ,  ,^1/2) 


2/3 


(1.3)' 


Using  the  Interaction  potential  (1.2),  Firsov  (Ref.  2)  obtains  an  expression  for 
which  he  approximates  by 


tr 


=  27re  M 


l“2  (\  ^2/':cm<“.  *  “2>)'  >'>  f  *  (I-'** 


where  is  the  energy  in  ev  in  the  center  of  mass  system: 


E _ =  ^E/M 


cm  ■  1 

^  =  Mj  M2/(M^  +  Mg) 


The  resultant  stopping  power'^  is 


=  1.3  X  10 


-13 


__LI2_ 

Mg  E/M 


-  (i . 


0.23 


M, 


Zj 


■) 


2 

cm  ev 


(1.5) 


To  obtain  the  energy  loss  due  to  excitation  and  ionization,  Firsov  (Ref.  3)  assumes 
that  there  is  a  friction-like  force  between  the  orbital  electrons  of  the  two  atoms  while 
they  are  passing  each  other.  The  resultant  conversion  of  the  kinetic  energy  of  rela¬ 
tive  motion  into  excitation  energy,  is  due  to  the  transfer  of  momentum  by  electrons 


^Stopping  power  is  the  average  energy  loss  per  incident  particle  per  centimeter  path 
in  a  gas  of  unit  number  density. 
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from  one  particle  to  the  other  in  the  region  where  the  electronic  shells  overlap.  The 
energy  lost  to  electron  excitation  or  ionization  at  a  given  impact  parameter  (cm) 
and  velocity  V  (cm/ sec)  is 

(Z,  +  X  4.3  X  lO'®  V 

€  =  - - - - ^ - ;r  (1.6) 

[l  +  3.1(Zj  +  10^  rJ 

Integrating  over  all  impact  parameters  to  obtain  the  total  stopping  power  gives 


Wg  =  2ir /  €  R^dR^  =  0.234  x  lO"^^  (Zj  +  Zglv  cm^  ev  (1.7) 


The  ionization  cross  section  at  a  velocity  v  is  given  by 


0-  =  (T, 


(1.8) 


with 


32.7  ,„-16  2 

- 273 


(Zi  +  Zg) 


and 


(1.9) 


23.3  £ 


^  = 


i  6 

•  'g^g  10  cm/  sec 
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in  (ev)  being  the  first  ionization  potential  of  the  gas  atom.  The  variation  of  a 
with  velocity  was  computed  from  Eq.  (1.8)  for  =  45  =  7: 


V  (cm/sec) 

2 

o(cm  X  10 

3  X 

10 

0.499 

5  X 

10® 

0.897 

8  X 

10® 

1.421 

1  X 

10^ 

1.732 

3  X 

lo”^ 

4.064 

5  X 

lO*^ 

5.748 

A  comparison  of  the  predictions  of  Eq.  (1.8)  with  experiment  has  been  carried  out  by 

0 

Fedorenko  (Ref.  4).  The  velocity  range  of  the  colliding  ions  varied  from  7  x  10  to 

7 

9  X  10  cm/sec,  and  in  this  region  there  is  agreement  between  theory  and  experiment 
to  within  a  factor  of  2.  Equation  (1.8)  is  even  useful  for  the  lighter  atoms, where  one 
finds  it  to  agree  to  within  a  factor  of  3.  A  comparison  of  the  theory  with  the  experi¬ 
mental  data  of  Ref.  5  for  electron  loss  cross  sections  for  lithium  in  helium: 

Li  +  He  — Li^  +  He  +  e"  ^  AK  (1.10) 

is  made  below 


(cm/ sec ) 

cr(cm^  X  10  ^®)  (Firsov) 

,  2  ,„-16.  , 
CT(cm  X  10  )  (expt) 

1  X  lo"^ 

1.45 

0.6 

6  X  10^ 

2.52 

1.0 

8  X  10^ 

3.54 

1.3 

A  comparison  of  the  stopping  powers  due  to  elastic  collisions  with  those  due  to  ioniz¬ 
ing  collisions  may  be  made  with  the  help  of  Eqs.  (1.5)  and  (1.7).  As  a  typical  example, 
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consider  a  fission  ion  of  atomic  mass  100  and  atomic  number  45.  The  stopping 
powers  are: 


v(cm/Bec) 
4  X  10^ 

1  X  10^ 


W g  ( cm^  ev )  (elastic) 

1.655  X  lO"^® 
1.79  X  lO"^® 


Wg  (cm®  ev )  (inelastic) 

4.69  X  10” 

-14 

1.17  X  10 


The  stopping  power  due  to  inelastic  scattering  falls  off  rapidly  with  decreasing  energy, 

7 

being  30  percent  of  that  due  to  elastic  scattering  at  a  velocity  of  4  x  10  cm/ sec,  and 

7 

only  7  percent  at  a  velocity  of  1  x  10  cm/sec. 


The  stopping  power  due  to  electron  capture  (charge  transfer)  is  more  difficult  to 
estimate.  For  lack  of  sufficient  theoretical  or  experimental  data,  it  is  necessary  to 
rely  on  empirical  rules.  From  an  examination  of  available  experimental  data  on 
charge  transfer  reactions, Hasted  (Ref.  6)  has  observed  that  the  maximum  charge- 
transfer  cross  section  occurs  at  a  veloci^  v^  of  the  incident  particle  given  by 


a  AE 
V  h 


m 


1 


(1.11) 


where  AE  is  the  energy  defect  measured  in  electron  volts  and  "a"  has  the  dimen¬ 
sion  of  length.  Hasted  claims  the  best  fit  is  given  by  a  =  8  A  • 

If  we  express  AE  in  ev  , 

v^  =  1.9  X  10^  AE  cm/ sec  (1-12) 

For  velocities  in  the  adiabatic  region,  the  charge  transfer  cross  section  varies 
according  to 


a  ~  c  exp(- k  a  AE/v) 


(1.13) 


where  c  and  k  are  constants  for  any  particular  reaction.  There  is,  of  course,  no 
rigorous  theoretical  basis  for  Eq.  (1.13).  For  velocities,  greater  than  that  for  which 
‘^max  ^  Increasing  energy  as  some  inverse  power  of  v  . 

For  a  typical  reaction 


D'*'  +  A^D  + 


the  energy  difference  will  be  of  the  order  of  7  ev  ,  so  that  the  velocity  at  which  a 


max 


occurs  is  approximately  1.3  x  10  cm/sec  .  It  would  thus  appear  that  is 

reached  for  velocities  well  outside  the  range  of  velocities  with  which  we  are  concerned. 


The  stopping  power  due  to  capture  is 


IEo 

^  n  c,n  Mj  ^  c,n 


For  a  velocity  of  4  x  10  cm/ sec  ,  E  =  83.6  kev 


and 


0. 455  ev 


_  2  7 

The  value  of  a  will  be  of  order  of  10  cm  and,  at  a  velocity  of  4  x  10  cm/ sec, 
max 

a  should  be  much  less  than  cr  according  to  Eq.  (1.13).  Nevertheless,  to  obtain 

mdx 

a  gross  overestimate  of  the  effect  of  charge  exchange  on  the  energy  loss,  let  us  assume 
7  “15  2 

that  at  V  =  4  X  10  cm/sec,  is  10  cm  .  Assuming  E^^  to  be  approximately 

10  ev  ,  we  find  the  stopping  power  due  to  charge  transfer  must  be  less  than 
14  2 

1  X  10  cm  ev  ,  which  is  less  than  that  due  to  inelastic  scattering.  Thus,  energy 
losses  due  to  charge  transfer  can  be  neglected. 
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Neglecting  the  energy  losses  due  to  ineUstic  collisions  altogether,  we  expand  the 
logarithm  in  Eq.  (1.4)  to  obtain 


a  1  . 

‘^tr  E 

X 


^tr 


1/v^ 


(1.14) 


From  Eqs.  (1.1)  and  (1.14) 


t  =  -KN  (1.15) 

where  K  is  a  function  of  A  and  Z  ,  and  N  is  the  number  density  of  the  atmosphere. 

Because  K  is  a  slowly  varying  function  of  its  arguments,  the  error  involved  in  work- 

-3 

ing  with  K  =  1.8  X  10  (the  value  for  a  typical  fragment  of  A  -  100  and  Z  =  45) 
will  not  be  great. 

9  with  the  vertical 
(1.16) 

Introducing  the  density  of  the  atmosphere  p  (z )  and 

N  =  ^  (1.17) 


For  a  particle  moving  along  a  straight  line  path  making  an  angle 

ds  =  ds  cos  6 


I  / 
I  / 
-T-V. 


where  A'  is  Avogadro*  s  number  and  M  the  mean  molecular  weight  of  the  atmos- 
idieric  atoms,  we  can  write  Eq.  (1.15) 


which  can  be  integrated  to  give 


Z 

0 

The  precise  form  of  the  stopping  law  used  depends  upon  the  assumptions  made 
concerning  p  (z )  in  evaluating  Eq.  (1. 19).  The  simplest,  and  the  one  adopted  here 
is  that  the  hydrostatic  condition 

dp  =  -  g  p(z)  dz  (1,20) 

is  satisfied.  Then,  if  we  neglect  the  variation  of  g  with  altitude,  Eq.  (1.19)  gives 
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Section  2 

INELASTIC  COLLISIONS  OF  FISSION  PARTICLES  IN  THE  ATMOSPHERE 

2. 1  THE  CHARGE  STATE  OF  A  BEAM  OF  PARTICLES 

According  to  the  results  of  the  preceding  section,  the  electron  capture  and  loss  cross 
secticms  do  not  contribute  to  the  energy  losses  of  particles  moving  throu^^  the  atmos- 
phere  at  velocities  equal  to  or  less  than  about  5  x  10  cm/sec.  Nonetheless,  these 
cross  sections  can  be  of  importance  in  determining  the  overall  charge  state  of  a  col¬ 
lection  of  particles. 

Let  and  n  denote  the  equilibrium  number  densities  of  the  debris  ions  and  atoms, 
respectively.  Neglecting  the  formation  of  more  highly  ionized  states  and  negative  icms 
of  the  debris  we  have  in  equilibrium 


“i'^c  ^ 

where  cr^  and  are  the  capture  and  loss  cross  sections.  The  charged  fraction  of 
the  beam  is 


f(i) 


n. 


n  +  n. 


1 


(2.2) 


For  a  typical  electron  loss  reaction 

D  +  A  —  +  e”  +  A  +  AE 

we  can  use  Eq.  (1. 8)  with  E^  of  the  order  of  7  ev.  The  capture  cress  section  can  be 
estimated  with  the  help  of  Hasted's  (Ref.  7)  empirical  relation  analogous  to  Eq.  (1. 13) 

=  c  exp  {-  aAE/4hv)  (2.3) 
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Taking  aAE/4h  =>;  3.3  x  10  for  a  typical  charge  transfer  reaction 

0“^  +  A  -  +  D 

we  have  the  following 

V( cm/sec)  o.  (cm^  x  10  a  (cm^  x  lo 

Mt  C 

4  X  10*^  7.64  .44 

1  X  10  2.96  .037 

We  see  that  a  falls  relatively  more  rapidly  with  energy  than  does  a  ,  and  over  the 

c  t 

range  of  velocities  of  interest  >  o’g  •  Hence,  from  Eq.  (2.2),  the  fraction  of  the 
beam  which  remains  charged  is  greater  than  1/2.  There  Is  some  experimental  con¬ 
firmation  of  this  result  (Ref.  8). 

2.2  RANDOM  INELASTIC  COLLISIONS 

in  addition  to  making  rough  estimates  like  Eq.  (2.  2)  possible,  the  inelastic  cross  sec¬ 
tions  determine  the  probability  for  a  particle  to  undergo  a  charge-changing  collision. 
To  determine  the  probability  distribution  for  such  a  collision,  we  consider  a  neutral 
particle  moving  through  a  medium  of  density  N  along  a  path  s  making  an  angle  0 
with  the  vertical.  Let  cr  be  the  cross-section  for  ionization  of  the  particle. 

The  change  in  the  number  n  of  neutral  particles  along  a  segment  of  path  ds  is  given 
by  the  well-known  attenuation  law 


dn  =  -  noNds 

and,  using  Eqs.  (1. 16)  and  (1. 17),  this  can  be  written 


dn 

n 


B 

cos  e 


p (z )  dz 


(2.4) 


(2.5) 
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with 


B 


_  cA* 
M 


(2.6) 


Integrating  Eq.  (2. 5)  from  an  arbitrary  Initial  altitude  to  an  arbitrary  final  altitude 
z  ,  we  have 


"/no  = 

where 


(2.7) 


z 

P(z,z^)  =  J  p(z)dz  (2.8) 

*o 

Because  COB  6  and  P(z,z^)  have  the  same  sign  and  are  negative  for  particles  moving 
downward,  n/n^  is  always  a  decreasing  exponential.  Observe  that  the  integral  Eq.  (2. 8) 
also  occurs  in  the  stopping  law,  Eq.  (1. 19). 

Let  n  be  the  number  of  particles  colliding  with  the  ambient  atoms  in  the  interval 
c 

( E  ,  z )  .  The  fraction 

'  o  •  • 


(2.9) 


is  the  probability  for  a  particle  to  have  an  icmizing  ooUislon  before  reaching  the  altitude 
z  .  Differentiating 


gives  the  spatial  rate  of  collision 


B 

cos  9 


exp 


1--^ 
I  cos  9 


P(*.*0>1  P(*)d* 


(2. 10) 


(2. 11) 
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The  density  p(z)  is  defined  in  0  <  z  <  <»  .  if  the  integrated  density  diverges,  thflt 
is  if 


OO 

P(«.Zo)  ^  /  P(z)dz  =  +  00  (2,12) 

z 

o 

then  Eq.  (2. 11)  is  a  normalized  probability  density  for  collisions  and  Eq.  (2. 10)  is  the 
corresponding  probability  distribution.  If  the  integrated  density  P(oo,z^)  is  finite, 
the  fraction  of  particles  escaping  to  infinity  without  undergoing  an  ionizing  collisicm  is 
finite; 


0  ' 

the  fraction  colliding  in  the  atmosphere  is  1  -  n  /n  ;  and  hence,  the  probability 

OO  Q 

distribution  for  collisions  in  the  open  interval  z^  <  z  <  «>  is 


1 


1  -  n 


(2. 14) 


The  standard  Monte  Carlo  procedure  is  to  choose  a  random  number  R  from  the  uniform 
distribution  on  (0,1)  and  to  solve  the  equation 


1  -  n/n 

"  '  1  -  .  /n  <2- 1=) 

00  0 

for  the  altitude  z  at  which  a  collision  occurs.  A  collection  of  altitudes  so  determined 
has  a  distribution  approximating  Eq.  (2. 14).  Calculating  z  directly  from  Eq.  (2. 15) 
can  lead  to  complicated  and  slow  computations  because  we  have  to  solve  the  equation 


/p(z)dz  =  In 


exp 


B 

cos  6 


f  P(z)dz 


(2.16) 


for  z  . 
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Alternatively,  we  can  interpret  Eq.  (2. 10)  aa  a  normalized  distribution  on  the  closed 
interval  z^  s  z  <  -t-  «  with  a  jump  discontinuity  at  the  point  z  =  +  <«  .  The  prob¬ 
ability  associated  with  z  =  +  <»  is  n  /n  .  With  the  discontinuous  distributicm,  the 
procedure  for  determining  z  is  to  solve 


for  z  if 


R  =  1  -  n/n 

o 


(2. 17) 


«  ^  1  <2. 18) 

holds  and  to  set  z  =-<-«>  if  Eq.  (2. 18)  fails.  Because  R  and  1  -  R  are  simul¬ 
taneously  uniformly  distributed  on  (0, 1).  we  solve,  in  practice, 

P(z.z^)  =  -  1„R  (2.19) 

for  z  if 


R  >  n^/n^  (2. 20) 

and  set  z  =  <«  if  Eq.  (2.20)  fails. 

For  a  finite  atmosphere,  one  in  which  p(z)  is  defined  for  0  :s  z  <  H,  we  use  the 
same  procedure, with  the  fraction  n„/n  replacing  the  fraction  n  /n  .  In  practice, 
we  solve  Eq.  (2. 19)  for  z  and  say  that  the  particle  has  reached  H  without  suffering 
a  collision  if  z  >  H  . 

2. 3  COLLISIONS  IN  VARIOUS  ATMOSPHERES 

The  explicit  calculation  of  the  collision  altitude  is  uniquely  determined  by  the  assumption 
made  on  the  behavior  of  p  ( z )  with  z  and  the  resultant  evaluation  of  Eq.  (2. 8). 


16 


2.3.1  Tabulated  Atmosphere 


Svqijpose  p(z)  is  known  as  a  table  of  values  of  p  vs.  z  in  an  interval  0  <  c  <  z  <  H  . 
By  a  suitable  interpolation  rule  we  can  determine  P(Zq)  for  any  z^  in  the  interval.  We 
can  also  construct  values  of  P(c,z)  and  P(H,z)  ,  and  by  interpolation,  those  of 
P(c,z^)  and  P(H,  z^)  .  Then  the  probability  of  escaping  over  the  iq;>per  boundary  H 
without  a  collision  between  the  starting  point  z^  and  H  ,  is 


V°o  ■  1-  (2.21) 

and  the  probability  distribution  for  collisions  in  the  interval  (H,  z^)  is 


1  -  "h^”o 


(2.  22) 


with  n/n^  given  by  Eq.  (2.7).  Equating  (2. 22)  to  a  random  number  R  and  solving 
for  z  ,  yields  a  collision  point  distributed  according  to  Eq.  (2. 14),  i.e. , 


z 

/  p  (z )  dz  =  -  In  [  1  -  R  ( 1  -  njj/n^)  ]  (2. 23) 

z 


Interpolation  among  the  tabulated  values  of  the  integrated  density  gives  z  . 

2. 3. 2  Constant  Atmosphere 

Let  p  ( z )  =  p  a  constant  in  the  Interval  0<c<z<H.  In  this  case 

P(z,z^)  =  p(z  -  z^)  (2.24) 

and  the  probability  of  escaping  over  the  upper  boundary  is 

njj/n^  =  exp  I  -  ^  p(H  -  z^)|  (2.25) 
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Because  the  integrated  density  diverges  at  infinity,  we  can  replace  a  finite  atnoosphere 
with  constant  density  an  infinite  atmosphere  with  the  same  constant  density.  In  such 
an  atmosphere,  the  probability  distribution  for  collisions  is  sinq>ly 


1  -  n/n 
o 

Choosing  a  random  number  R  and  solving  for  z,  we  get 


z 

If  this  value  of  z  fails  to  satisfy 


cos  6 
PB 


InR 


(2.26) 


(2. 27) 


c  <  z  <  H 

it  is  rejected  as  a  possible  collision  point.  Note  that  the  probability  of  rejecting  such 
a  z  is 


/  <“«>(■  -  Vl  |-  <«  -  Vl 

which  is  precisely  the  probability  of  escape  over  the  upper  boundary  H  . 

2. 3. 3  Linear  Atmosphere 

Suppose  the  density  decreases  linearly  according  to 

p(z)  =  b  -  az  ,  a  >  0  (2.29) 

lnO<c<z<H.  Because  p  ( z )  >  0  ,  we  have 

b/a  >  H  (2.30) 
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The  integrated  density  is 


P(z.z^)  =  (z  -  z^)  [b  -  I  (z  +  z^)] 


(2.31) 


a  quadratic  polynomial  with  zeros  at  and 


2b 

z  =  —  -  z 
a  o 


The  fact  that  z^  <  H  ,  and  Eq.  (2. 30)  imply 


^-z  >  H 
a  o 


(2. 32) 


(2.33) 


Hence,  the  second  zero  occurs  above  H.  The  collision  point  is  determined  in  the 
usual  manner. 

2.3.4  Exponentially  Decreasing  Atmosphere 


Suppose  the  density  is  given  by 

p ( z )  =  a  exp  I  -  b ( z  -  c )  | 

with  a,  b>0in0<c<z<H.  The  integrated  density  becomes 


(2.34) 


PCz.Zq) 


=  I  exp  |-  b(z^  -  c))  -  exp  |-  b(z  -  c))]  =  q(z^)  -  q(z) 


(2. 35) 


with 


q(z)  =1  exp  {-  b(z  -  c)) 

The  probability  of  escape  over  the  upper  boundary  H  is  now 

=  exp  |- |q(x^)  -  ,(H)|| 


(2. 36) 


19 


and  solving  the  equation 


1  -  n/n 

for  z  gives 

where 


p  =  _  in[i  .  R(i  _  n^/n^)] 


2. 3 .  S  Hydrostatic  Atmosphere 

The  previous  assumptions  concerning  p  ( z )  are  either  manifestly  unsuitable  or  lead 
to  computations  which  are  too  slow  for  Monte  Carlo  calculations.  The  hydrostatic  re¬ 
lation,  Eq.  (1.20) 


dp(z)  =  -  gp(z)  dz 

where  p  ( z )  is  the  pressure  at  altitude  z  and  g  is  the  acceleration  due  to  gravity 
can  be  Integrated  to  give 

P  -  P 

P(z,z^)  =  -  (2.38) 

where  g  is  assumed  constant  over  the  altitudes  of  interest.  The  fraction  surviving 
Eq.  (2.7)  becomes 


n/n^  =  exp 
0 


B 

cos  9 


(2.39) 
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Using  the  alternative  procedure,  Eqs.  (2. 19)  and  (2.20),  and  solving  for  the  collision 
pressure  rather  than  the  collision  altitude  gives 

P  =  Pq  +  g  In  R  (2. 40) 

as  a  provisional  value  of  p  to  be  accepted  only  if 

P  >  Pjj  (2. 41) 

Interpolation  in  the  table  of  p  vs  z  (Section  5)  determines  the  collision  altitude.  The 
altitudes  so  determined  have 


1  -  n/n 

as  a  normalized  probability  distribution. 

The  assumptions  usually  made  in  integrating  the  hydrostatic  relation  are  (i) 

P(z)  =  ^  P(z)  (2.43) 

where  k  is  Boltzmann's  constant,  M  the  molecular  weight,  and  T  the  temperature; 
(ii)  the  variation  of  g,  M,  and  T  can  be  neglected.  In  the  presence  of  these  assump¬ 
tions,  we  easily  find  that 

P(z)  =  P(Zo)  I"  ■  *0^1 

an  exponentially  decreasing  atmosphere,  |  cf  Eq.  (2. 34)  ] . 
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2.4  COLLISIONS  ALONG  A  SPIRAL  PATH 


Let  the  z'  axis  of  a  new  rectangular  coordinate  system  have  a  polar  an^  8^  with 
respect  to  the  z  axis.  Suppose  a  particle  travels  akmg  a  spiral  path,  say  the  helix 


X*  =  r  cos  u 

y»  =  r  sin  u  (2. 45) 

z*  =  bu 


where  r  is  the  radius  of  gyration  and  2xb  is  the  distance  the  particle  advances  almg 
the  z*  axis  in  aae  revolution.  If  we  unroll  the  path  for  a  single  revolution,  we  get  the 
following  sketch 


Here  6'  is  the  pitch  angle  of  the  spiral.  We  see  that  the  parameter  u  is  connected 
with  the  arc  length  s  by  the  relation 


u  =  s/(r^  +  b^) 

and  that 

/  2  2\^/2 
cos  O'  =  b/^r^  °  ) 

holds.  Thus,  we  have 


(2.46) 


(2.47) 


or 


1/2 


=  bs/^r^  *  ® 


(2.48) 


dz*  »  cos  0'  ds 


(2.49) 
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Moreover,  because  the  polar  angle  of  the  z'  axis  is  6^  we  have,  with  respect  to  the 
z  axis, 


dz  =  cos  0^  dz*  (2.50) 

Combining  Eqs.  (2.49)  and  (2.50)  gives 

dz  =  cos  6^  cos  0*  ds  (2.51) 

Comparing  Eq.  (2. 51)  with  Eq.  (1. 16)  shows  that  the  analogue  for  Eq.  (2. 5)  in  the 
case  of  a  spiral  path  is 


dn 

n 


B 

cos  cos  6' 


p(z)  dz 


(2. 52) 


Thus,  we  may  simply  replace  cos  8  by  the  product 

cos  9^  cos  0'  (2. 53) 

in  the  formulas  for  straight  line  paths  to  obtain  the  formulas  appropriate  for  spiral 
paths. 


We  also  make  the  same  replacement  when  applying  Eq,  (1. 21)  to  describe  the  slowing 
down  along  a  spiral  path. 
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SectloD  3 

BETA  DECAY  OF  FISSION  FRAGMENTS 


A  correct  statistical  treatment  of  beta  decay  would  require  a  knowledge  of  all  decay 
chains  and  the  lifetimes  and  branching  ratios  of  all  the  elements  in  these  chains. 

Because  such  detailed  data  are  not  known  e^eri mentally,  we  are  forced  to  make  model 
assumptions,  hi  this  we  are  guided  by  a  few  well-estaUished  experimental  facts  about 
beta  decay.  According  to  Way  and  Wigner  (Ref.  9)  the  combined  beta  activity  of  all  fis¬ 
sion  products  and  succeeding  generations  stays  roughly  constant  for  1  sec  and  then  drops 
_1  2 

as  t  '  .  The  length  of  the  decay  chains  varies,  but  on  the  average  about  three  betas 
are  emitted  per  fission  product.  We  define  the  beta  activity  A(t)  to  be  the  average  num¬ 
ber  of  betas  per  fission  product  per  unit  time  so  that 

00 

/A(t)dt  =  3  (3.1) 


With  this  normalization,  the  Way- Wigner  law  takes  the  form 


A(t) 


its 
2’  * 

1.-1. 2 

2 


1  sec 

,  t  >  1  sec 


(3.2) 


In  our  model  we  assume  that  all  fission  products  produce  exactly  three  electrons.  For 
the  Monte  Carlo  method,  we  need  probability  distributions  for  each  of  the  three  decays. 

Let  us  Introduce  the  fractions  of  fission  products  n^,  n^,  and  n^,  which  are  present 
after  0,  1,  2,  and  3  decays  have  taken  place.  At  t  =  0,  no  decays  have  occurred,  so 
that  the  initial  conditions  are 

n^  =  1  and  n^^  =  Ug  =  ng  =  0  (3.3) 
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Let  us  introduce  three  positive  functions,  fj^(t) ,  f2(t) »  f3(f ) .  describing  the 
rate  of  production  of  betas  by  first,  second,  and  third  decays.  The  fractions  n^  obey 
differential  equations  of  the  form 

dn 

-df  ■  - 


dt 


and  the  normalization  condition  for  all  times 


(3.4) 


=  1  -  n. 


"l  ■  "2 


(3.5) 


Because  there  are  not  enough  experimental  data  to  determine  the  three  rate  func¬ 
tions,  we  can  make  the  assumption  that  they  have  a  common  value,  say  f(t). 


Let  us  introduce  the  scaled  time 


t 

s  =  J  f(u)du 
o 

The  system  [Eq.  (3.4)]simpllfies  to 


%  _  „ 

ds  0 


"0  ■  “1 


(3.6) 


(3.7) 
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which  has  the  solution 


satisfying  the  initial  conditions  required  by  Eq.  (3.3). 


We  note  that 


J  nj  ( a )  ds  =  1  , 

o 


i  =  0,  1,  2, 


(3.8) 


(3.9) 


holds  for  each  of  Eqs.  (3. 8).  In  terms  of  the  scaled  time  s  ,  the  analogue  of  the  rate 
function  f  is  unity.  Consequently,  we  can  identify  n^  ,  i  =  0,  1,  2  as  the  probability 
density  of  the  ( i  +  1  )st  decay  time  In  scaled  time.  The  corresponding  densities  in  real 
time  are 


t 

f(t)nj  j  f(u)du  ,  i  =  0,  1,  2  (3.10) 

o 

Because  first,  second,  and  third  decays  for  an  individual  particle  occur  sequentially, 
we  need  a  procedure  for  choosing  three  random  vairables  >  ^2  ’  *"3 

0  <  ti  <  tg  <  t3  (3. 11) 

such  that  t^  is  a  sample  from  the  distribution  of  i^  decay  times  Eq.  3. 10  .  Out 
strategy  is  to  make  use  of  the  one--one,  monotone,  and  continuous  relation  between  s  and 
t,  Eq.  (3.6),  to  solve  the  analogous  problem  in  scaled  time,  and  then  to  transform  the 
random  variable  s^  into  random  variables  t^. 
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If  ,  i  =  1,  2,  3,  are  independent  random  variables  each  with  e  as  probability 

density,  then  the  random  variable  s,,  =  u-  +  u„  has  n,  (s)  =  se  ®  as  probability 

^  ^  ^  ^  /  2  \  ■ 
density,  and  the  random  vairable  Sg  =  Uj  +  Ug  +  =  Sg  +  Ug  has  Og  (s)  =  (s  /2je 

as  probability  density.  Consequently,  we  solve  the  problem  in  scaled  time  by  choosing 

three  random  numbers,  ,  i  =  1 , 2 , 3  ,  from  the  uniform  distribution  on  (0,1),  by 

setting  Uj  =  -InR^ ,  and  by  taking 

Sj^  =  Uj  =  -  InRj 

®2  "  “2  "  ®1  ■  ^"^2  (3.12) 

S3  =  “1  ^  “2  ^  “3  =  ®2  ■  ^3 


Requiring  the  total  activity  resulting  from  all  three  generations  to  reproduce  the  Way- 
Wigner  law  yields 


A  (t)  =  f  (t)  -i  -f  ng  j  = 


ds 

dt 


l  +  s+^ 


(3. 13) 


Integrating  Eq.  (3.  13)  with  respect  to  t  gives 


t 

j  A  (t)  dt 
o 


2 

1  +  s  +  Y 


-  s 


e 


ds 


(3. 14) 


Under  the  hypothesis  that  the  Way-Wigner  law  holds,  the  relation  Eq.  (3. 14)  between 
s  and  t  is  equivalent  to  Eq.  (3. 6)  but  does  not  involve  a  knowledge  of  f  ( t ) .  Carrying 
out  the  integrations  in  Eq.  (3. 14)  and  solving  for  t  gives 
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28 


0  s  B  <  1/2 


with 


1/2  «•  S 


S  =  3  - 


^3  +  28 


-B 


(3. 15) 


Combining  Eqb.  (3. 12)  and  (3. 15)  given  decay  times  t^  ,  i  =  1,  2,  3,  satiefying  (3. 11), 

The  uniform  rate  function  f  ( t )  is  nearly  conatant  for  0  s  t  ^  1  and  decreaaes  for 
t  1  .  In  fact,  manipulating  Eq.  (3. 14)  yields 


2^ 


2+28+8 


f(t) 


1  + 


4  »  88 

a  +  I8  +  S 


it 


t  <  1  860, 


t  >  1  sec 


with  the  value  of  s  for  which  t  =  1  being  approximately  0.51. 


(3. 16) 


A  simpler,  alternative  way  to  obtain  three  decay  times  for  a  particle  is  to  choose  three 
samples,  t  ,  from  the  distribution  with  the  Wt^-Wlgner  law  as  probability  density,  to 
order  the  samples  so  that  Eq.  (3. 11)  holds,  and  then  to  declare  the  resulting  tj  to  be 
the  1  decay  time.  This  procedure  amounts  to  having  different  rate  functions  in  the 
differential  |Eq.  (3.4))  . 
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Interpreted  as  a  probability  density,  the  Way-Wigner  law  takes  the  form 


w(t)  = 


1/6  ,  t  1  sec 


.  t  >  1  sec 


(3. 17) 


and  the  corresponding  probability  distribution  is 


W(t) 


t/6  ,  t  1  sec 

1  -  5t“®‘  ^/6  ,  t  >  1  sec 


(3.18) 


To  obtain  a  sample  from  this  distribution,  we  choose  a  random  number,  R  ,  from 
the  uniform  distribution  on  ( 0, 1 )  and  set 


t  =  W’^R) 


(3.19) 


Because  R  =  W(t)  is  a  monotone,  increasing,  and  coatinuous  function,  ordering 
the  three  random  numbers,  R  ,  to  satisfy 

0  <  Rj  <  Rg  <  Rg  (3.20) 

is  equivalent  to  ordering  the  three  corresponding  samples,  t  ,  to  satisfy  Eq.  (3. 11). 
The  probability  densities  for  the  minimum,  middle,  and  maximiun  of  three  independent 
random  numbers  are 


p^(R)  =  3(1  -  R)^ 

P2(R)  =  6R(1  -  R) 

Pg(R)  =  3R^  (3-21) 
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We  regard  S  *  W(t)  as  a  mailing  from  real  time,  t  ,  to  another  scale  time,  R  . 
In  the  alternative  procedure,  the  Pj(R)  ,  i  1, 2, 3  ,  are  the  densities  of  i^  decays 

in  scale  time,  R  .  The  corresponding  rates  at  which  the  fractions  n^  ,  n^  ,  and 
n^  decrease  are 


dno 

dn^ 

=  p^(R)  -  P2(R)  (3.22) 

dn„ 

=  PjCR)  -  P3(R) 

Integrating  Eq.  (3. 22)  with  the  initial  conditions  (3.3)  gives 

Dq  »  (1  -  R)® 

nj  =  3R(1  -  R)^  (3.23) 

Ug  =  3R*(1  -  R) 

^y  introducing  three  different  rate  functions 

*2*®*  "  1  -  R 

g3(R)  =  (3.24) 
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we  can  rewrite  Eq.  (3.  22)  as 


dn 

=  -  «!<«> 
dn^ 

dif  ~  "o  “  ^2^^^  “l 

dn« 

-g^-  =  g2(R)  “i  -  KsCR)  ^2 


which  is  an  analog  of  the  systemfEq.  (3.7)]in  the  first  model. 

By  using  Eq.  (3. 19)  to  transform  Eq.  (3. 25)  from  scale  time,  R  ,  to  real  time,  t  , 
we  get  a  system  of  differential  equations  of  the  form  Eq.  (3.4)  in  which  the  rate 
functions  are 

V‘>  = 

^3(^)  ~  0  _  ^  (3.26) 

for  0  ^  t  s  1  and 

f^(t)  =  3/5t 
f2(t)  =  2/5t 

fg(t)  =  l/5t  •  (3.27) 

for  t  >  1  .  These  rate  functions  increase  for  0  :£  t  ^  1  and  decrease  as  1/t  for 
t>  1  . 
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We  adopt  the  procedure  of  Eqs.  (3. 18,  19,  and  20)a8  the  model  for  beta  decays  pri¬ 
marily  because  it  leads  to  simpler  computations.  Moreover,  the  rate  of  successive 
decays  decreases.  For  a  typical  fission  fragment  this  behavior  of  decay  rates  seems 
intuitively  plausible  because  of  the  excess  of  the  neutron-to-proton  ratio  of  the  particle, 
which  gives  rise  to  its  instability,  decreases  with  successive  decays. 


The  graphs  of  injection  density  in  Ref.  10  are  based  upon  still  another  prescription  for 
the  decay  times.  In  those  computations  we  used  the  relation 


.. 

5 

5  exp  s 

s  s 

2 

.6  +  4s  +  s 

*  U; 

.0.4  +  lOs^  1  +  lOs^. 

(3.28) 


with 


8  =  Sj  =  -  In  Rj  ,  i  =  1,  2,  3  (3.29) 

to  compute  the  decay  times.  Ordering  the  random  numbers  Rj  ,  chosen  from  the  uni¬ 
form  distribution  on  (0,  1),  to  satisfy 

0  <  Rg  <  Rg  <  Rj  <  1  (3.30) 

insures  that  the  t.  increase  with  i  .  Equation  (3.  28)  is  an  analytic  approximation  for 
the  two  relations  in  Eq.  (3. 15). 

For  a  triple  of  random  numbers  ordered  as  in  Eq.  (3.  30),  s„  =  -  In  R_  is  the  largest 

v  U 

value  produced  by  Eq.  (3.29);  for  the  same  triple,  Sg  =  -  (In  R^  +  In  Rg  +  In  Rg)  is 
the  largest  value  produced  by  Eq.  (3. 12).  Thus  it  is  clear  that  the  third  decay  times 
produced  by  the  prescription  of  Eqs.  (3. 28),  (3. 29),  and  (3.  30)  must  be  smaller  than 
those  produced  by  the  prescription  of  Eqs.  (3. 12)  and  (3. 15). 
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What  is  true  for  third  decays  holds  in  general:  the  total  distribution  of  decay  times, 
without  regard  for  decay  number,  produced  by  the  present  prescription  has  many  more 
early  times  and  for  fewer  late  times  than  do  the  distributions  produced  by  the  two  pre¬ 
vious  prescriptions,  which  approximate  the  Way-Wigner  law. 

Figures  1  and  2  show,  for  two  sets  of  parameters,  the  variation  of  injection  density 
with  distance  for  the  three  prescriptions  for  the  decay  times;  Figures  3  and  4  show 

4 

the  variation  of  injection  density  with  altitude  for  the  same  data.  With  2. 5  x  10 
histories,  the  results  obtained  from  Eqs.  (3. 12)  and  (3. 15)  are  almost  identical  with 
those  obtained  from  Eqs.  (3. 18),  (3. 19),  and  (3.  20).  The  effect  of  the  prependerance 
of  early  decay  times  in  the  third  prescription,  Eqs.  (3.28),  (3.29),  and  (3.30),  is  to 
suppress  the  tails  and  to  accentuate  the  central  peak  in  the  distributions  of  injection 
density. 

We  observe  that  the  parameter  change  (increasing  the  velocity)  broadens  the  central 
peak  for  all  prescriptions  and  that  the  broadening  is  greatest  for  the  third  prescription. 

The  discussion  of  the  relative  influence  of  the  various  parameters  in  Section  3  of  Ref.  1 0 
is  based  upon  the  comparative  shapes  of  the  injection  density  curves  and  is  facilitated  by 
using  the  third  prescription  to  magnify  these  effects.  On  the  other  hand  effects  which 
occur  at  late  times  are  minimized  by  the  third  prescription.  Our  computations  suggest 
that  the  assumption  of  a  uniform  decay  rate,  the  first  prescription,  also  tends  to  obscure 
effects  occurring  at  late  times.  Consequently,  for  statistics  concerning  effects  at  late 
times,  we  recommend  the  second  prescription  for  decay  times,  Eqs.  (3. 18).  (3. 19),  and 
(3.  20).  This  second  prescription  is  in  the  current  version  of  the  code. 


33 


Eqs.  (3.12)  to  (3.15) 
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Fig.  1  Variation  of  Injection  Density  With  Distance  for  Different  Prescriptions  of  Decay  Times 
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Fig.  2  Vaxiati(Mi  of  Injection  Density  With  Distance  for  Different  Prescriptions  of  Decay  Times 
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for  Different  Prescriptions  of  Decay  Times 


32001-  \l 


(SilNH  AUVdllSdV )  *u 
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Fig.  4  Variation  of  Injection  Density  With  Altitude  for  Different  Prescriptions  of  Decay  Times 


Section  4 

IBM  7090  CODE  FOR  THE  MONTE  CARLO  MODEL 


The  fraction  of  the  beam  which  is  ionized  is  a  free  parameter  and  constitutes  one  of 
the  input  constants  of  the  problem.  A  random  number  chosen  between  0  and  1, 
compared  with  this  fraction,  decides  whether  the  particle  considered  is  neutral  or 
ionized.  From  then  on,  a  special  counter  (I)  keeps  track  of  the  charge  on  the 
particle. 

The  three  decay  times  are  computed  from  three  random  numbers  chosen  independ¬ 
ently  from  a  uniform  distribution  on  (0, 1)  and  then  ordered  to  satisfy 

0  <  R^  <  Rg  <  Rg  <  1  (4.1) 

The  decay  times  are  obtained  by  the  prescription,  derived  from  Eqs.  (3. 18) 
and  (3. 17), 


t 


i 


[6(1  -Rj)J 


if  Rj^s  1/6 
if  Rj>  1/6 


(4.2) 
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A  particle  initially  neutral  is  taken  as  traveling  in  a  straight  line  in  a  direction 
chosen  from  a  uniform  distribution  on  a  hemisphere.  Figure  1  shows  the  coordinates 
of  the  particle  with  9  and  <p  taken  as  uniformly  distributed  between  0  and  ir  ;  and 
0  and  2n  ,  respectively.  On  ionization  by  collision  or  beta  decay,  the  particle  is 
constrained  to  spiral  about  the  magnetic  field  line  passing  through  the  point  at  which 
ionization  occurs,  thus  neglecting  the  lateral  displacement  equal  to  the  radius  of 
gyration.  Because  the  latter  is  an  order  of  magnitude  lower  than  the  mean  free  path, 
the  resultant  error  in  the  coordinates  of  the  fragment  will,  on  the  average,  be  neg¬ 
ligible.  To  follow  the  motion  of  the  ionized  particles,  we  go  to  a  coordinate  system 
in  which  the  z'-axis  is  parallel  to  the  magnetic  field  through  the  origin.  In  this 
(magnetic  coordinate)  system,  the  polar  angle  of  the  particle  position  d'  is  also  its 
pitch  angle.*  For  the  geometry  of  Fig.  5  we  get 

cos  o'  =  cos  9  cos  9  -  sin  0  sin  9  sin  a  (4. 3) 

cos  (p  =  (sin  d/sin  0')  cos  (p  (4-4) 

Also,  from  the  figure,  we  see  that  if  the  particle  moves  a  distance  ds  along  the 
spiral  it  will  rise  through  a  verticle  dist^mce  given  by 

dz  =  ds  cos  0' cos  0  (4.5) 

This  equation,  which  is  the  same  as  Eq.  (2.51),  is  used  to  obtain  the  effective  path 
length,  because,  in  spiralling,  the  ionized  particle  has  a  greater  probability  of 
collision  than  a  neutral  particle  moving  through  the  same  verticle  distance.  The 
ionized  particle  remains  on  the  spiral  until  it  is  slowed  down  to  rest,  neutralized 
by  collision,  or  ceases  to  be  of  interest  because  it  has  decayed  three  times.  If  it 
remains  on  the  spiral,  it  is  followed  to  the  conjugate  point.  The  coordinates  of 
every  decay  point  are  noted,  but  the  motion  of  the  ionized  particle  remains 
unaffected. 


♦For  initially  ionized  particles,  0'  is  chosen  from  a  imiform  distribution  on  a  hemi¬ 
sphere  with  the  polar  axis  parallel  to  the  magnetic  field,and  the  calculation  proceeds 
as  above. 
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Throughout  its  passage  through  the  atmosphere  the  particle's  velocity  is  reduced 

according  to  Eq.  (1. 21)  until  it  escapes  with  the  velocity  it  has  at  the  altitude  H.  Let 

V  , ,  V  and  p  . ,  p  be  the  values  of  the  velocity  and  pressure  at  any  two  con- 
n-in  n-in  j 

secutive  points  along  the  path.  Generalizing  Eq.  (1.21),  we  have 


v 


2 


n 


2  2kA' 

Y  ^  . . 

n  -  1  Mg  cos  0 


-  Pn> 


Similarly,  generalizing  Eq.  (2.40)  yields 


(4.5) 


.  Mg  C0..5 
n  '^n-l  CTjA' 


(4.6) 


where  use  has  been  made  of  Eq.  (2.6).  Combining  Eqs.  (4.5)  and  (4.6)  gives  the 
alternative  form  for  the  slowing  down  law: 


2  2  ^  2k  ,  _ 
V.,  =  —  In  R 

n  n-1 


In  the  code  we  use  the  notation. 


(4.7) 


Mg  cos  0  _  * 
A'a.  i  ’ 


2k  A' 
Mg  cos  0 


=  B 


A  2k 
and  —  = 

•^i 


(4.8) 


We  note  the  appearance  of  the  inelastic  collision  cross  section  in  Eq.  (4.  7), 
although  the  corresponding  collisions  are  ignored  as  sources  of  energy  loss.  The 
explanation  is  immediate  if  we  remember  that  the  slowing  down  of  a  particle  must  be 
proportional  to  its  path  length.  The  effect  of  different  tjrpes  of  inelastic  collision  cross 
sections  can  be  studied  by  varying  the  constants  Aj  and  .  The  effect  of  various 
elastic  collision  cross  sections  could  be  investigated  by  varying  k  in  the  constants 
B  and  Cj  . 


A  collision  with  an  atmospheric  atom  reduces  the  charge  on  the  particle  by  unity,  the 
probability  for  electron  capture  being  definitely  greater  than  that  for  loss  because  of 
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the  relatively  high  valuee  of  the  ionlzaticm  potentials  beyond  the  first.  The  motion  of  the 
neutralized  particle  la  followed  by  reverting  to  the  unprimed  coordinate  system  by  means 
of  the  transformation  inverse  to  Eqs.  (4. 3)  and  (4.4),  viz. , 

cos  0  =  cos  O'  cos  0  -t-  sin  d'  sin  6  sin  O'  ,  (4.9) 

0  0 

cos  9  =  (sin  9 '/sin  9)  cos  9*  (4.10) 

where  9'  is  the  azimuthal  angle  in  the  magnetic  coordinate  system  and  is  chosen  from 
a  uniform  diatribution  between  0  and  r  . 

Decays  and  Inelastic  collisions  are  Independent, and  therefore  the  time  intervals  to 
decays  and  collisions  are  used  to  determine  which  event  occurs  first.  By  means  of 
a  random  number,  the  pressure  at  the  point  of  collision  is  determined;  the  pressure 
is  a  known  function  of  the  altitude  which  can  be  determined  from  the  tables  reproduced 
for  convenience,  in  Secticm  6.  Hence,  the  distance  travelled  and  the  time  to  cover  this 
distance  are  found,  and  this  time  is  compared  with  the  time  to  first  decay.  If  the  time 
to  decay  is  shorter,  the  time  to  collision  is  reduced  by  this  amount  and  the  residue 
compared  with  the  time  to  seccmd  decay.  If  the  time  to  collision  is  shorter,  the  time 
to  first  decay  is  reduced,  and  a  new  collision  time  is  calculated  and  compared  with  it. 

The  need  for  two  slowing  down  equations  arises  aa  follows.  Consider  a  neutral  particle 
at  a  starting  point  ^  and  suppose  it  is  to  undergo  an  ionizing  collision  at  a  point 
z^  .  We  determine  p^  and  hence  v^  by  either  Eq.  (4. 5)orEq.  (4.7).  Stqipose,  how¬ 
ever,  we  find  that  a  decay  will  occ\ir  enroute  to  z  ,  at  zl  .  There  is  no  direct 

n  n 

method  of  determining  this  point  and  so  we  proceed  in  the  following  way:  we  calculate 
the  average  velocity  from  2^  _  ^  to  z^  and,  with  the  time  to  decay,  find  the  path  length 
to  z^  and  hence  z^  and  pj^ .  To  proceed  with  the  history  we  require  v^  .  It  would 
clearly  be  very  inaccurate  to  use  v^  or  the  above  average  because  the  point  z^  can  be 
far  below  z^  and  very  close  to  2^  _  ^  .  Equation  (4. 7)  is  now  useless  and  we  have  to 
use  Eq.  (4. 6). 
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At  each  calculation  of  the  collision  time,  the  new  pressure  is  compared  with  that  at 
the  altitude  H  to  determine  whether  the  particle  escapes  from  the  atmosphere.  Neu¬ 
tral  particles  which  would  escape  are  checked  for  possible  decays  before  they  could 
do  so,  because  in  case  of  decay,  their  path  in  the  atmosphere  is  lengthened,  they  are 
slowed  down  more,  and  might  even  be  stopped.  Neutral  particles  which  do  escape  are 
checked  for  decays  up  to  the  altitude  of  20,000  km  above  which  they  are  considered 
to  be  lost. 

The  calculation  thus  proceeds  step  by  step  until  the  history  of  a  particle  is  terminated 
by  one  of  the  following: 

•  three  decays 

•  loss  at  the  conjugate  point 

•  loss  in  atmosphere  by  slowing  down  to  rest 

•  loss  in  atmosphere  by  scattering  in  the  backward  direction 
(cos  i9  or  cos  9'  negative  in  Eq.  (4.4)  or  (4.9) 

•  loss  by  escaping  from  the  atmosphere  and  traveling  out  of  region 
of  interest  (neutral  or  neutralized  particles  only) 

The  code  is  set  out  in  eight  blocks  of  instructions,  each  designed  to  cover  a  particular 
aspect  of  a  particle's  history  according  to  the  following  scheme: 


Instruction  Set  No. 

Function 

100 

Initialization 

200  Collision  Routine  for  Neutral  Particles  (CN) 

300  Collision  Routine  for  Ionized  Particles  (Cl) 

400  Decay  Routine  for  Neutral  Particles  (DN) 

500  Decay  Routine  for  Ionized  Particels  ( DI) 

600  Escape  Routine  for  Neutral  Particles  (EN) 

800  Escape  Routine  for  Ionized  Particles  (El) 

6600  Escape  Routine  for  Neutral  Particles 

with  zero  ionization  cross  section  ( EN° ) 


The  details  of  all  the  steps  involved  are  shown  in  the  flow  charts  of  Figs.  6  through  16. 
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BCAK  Of  NfijTtAU  MTH 

ENTIi^  DECAY  or  lOMZfD  DfCAY  OP  NIUTtAU  2«0  tOMlZAHON  CtOSS-SKTlON  ESCAff  OP  NiUTtACS  ESCA«  OP  IONIZED 
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FIk.  7  Master  Flow  Chart 


illilii 


START  OF  PROILEM  SET 


START  OF  A  NEW  HISTORY 


Fig.  8  biltlalization  Flow  Chart 


CALCULATION  OF  DECAY  TIMES 
Eqs.  (3.28,  29,  and  30) 


Fig.  9a  Decay  Times  Flow  Chart  (used  in  production  runs  for  Ref,  10) 
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CALCULATION  OP  DECAY  TIMES 
<Eq«.  (3. 18,  19.  and  30) 


Fig.  9b  Decay  Times  Flow  Chart  Yielding  Correct  Late  Time  Effects 
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CALCULATION  OF  REDUCED  TIMES  AND 
INITIAL  CONDITIONS 


117 


Fig.  9c  Decay  Times  Flow  Chart  (completion) 
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VN  =  SQRTF  (T) 

VB  =  0.5*(VN  +  VNB) 
RTC  =  DSO/VB 


Fig.  10  Neutral  Collision  Flow  Chart 
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Fig.  11  Ionized  Collision  Flow  Chart 
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CALCUWTION  Of  decays  OF  NEUTRAL  PARTICLES  I 


RTO  -  RTC 


XNI >  XN 
YNR  -  YN 
ZN»  -  4N 
PNB  •  PN 

VNE  •  VN 
K  •  K  I 
I  >  I  -I-  I 

RTD  =  RTD  -  RTC 
CTETAP  .  CTETA  •  CTETAO 
-STETA-SPHI-STETAC 


NE-NErl 
L=  L»  1 
I  •  I  ♦  I 

OSO  •  RTD  *  V8 
DZO  -  DSO  •  CTETA 
DXO  =  DSO  *  STETA  •  CPHI 
DYO  *  DSO  ♦  STETA  •  SPHI 
XN  XNS*  DXO 
YN  -  YNB  *  DYO 
ZN=  i.NB*OZO 


INTERPOLATE  FOR  PRESSURE 
IN  ALTITUDE  TABLE  iSSOO) 


5TETAP  =  5QRTF  (1 .0  -  CTETAP  *•  2) 
RTC  •  RTC  -  RTD 


RTD  =  TDTWO  \ 


Jrtd^tdthr 


Fig.  12  Neutrsd  Decay  Flow  Chart 
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Fig.  14  Neutral  Escape  Flow  Chart 
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Fig.  15  Ionized  Escape  Flow  Chart 
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ESCAPE  ROUTINE  FOR  NEUTRALS  WITH  ZERO  IONISATION  CROSS-SECTION 


Fig.  16  Neutral  Escape  (<T|  =0)  Flowchart 
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Section  5 

FORTRAN  LISTING  AND  SAMPLE  PROBLEM 

This  is  the  final  version  of  the  code  using  Eqs.  (3. 18),  (3. 19),  and  (ji.  20)  for  the 
decay  times.  This  version  gives  correct  late  time  effects. 


FORTRAN 

207-89  ORMONDE  X49343 
CHANGED  PROGRAM  TO  USE  TWO 
TIVE  MATRIX  NONE  ZERO 
DEMO  WITH  both  plot  AND 

NEw’^nSmE  ^Stb^NAMl®* 


KCT 

JUMP 

NO 

r.i 

(CCS 

CTO 

STO 

CTP 

STP 

CT 

ST 

CP 

SP 

CPP 

SPP 


(COUNTY 


STETAO 

CTETAP 

STETAP 

CTETA 

STETA 

CPHI 

SPHj 

CPHIP 

SPHIP 


DEMO  FINAL  VERSION 

SIGMAS  FOR  EACH  CASE  W/COLUMN  ONE  OF  NECA- 
AUGUST  14*  1962 

MATRIX  S-31-62  BACKWARD  TEST  IS  IN 

OES^AyPTliiJ®**^*’  ORMONDE  92-10  X45219 

total  CARD  COUNT 

INDICATOR-CROSS-SECTION 

NUMBER  OF  neutral  PARTICLES  IN  ZERO) 

?^IE"«85efS'SI'lA«iEtft' 

|isiNE^?H€fS%RIME 

sTn  theta 

COSIN  PHI 
SIN  PHI 

COSINE  PHI  PRIME 
SIN  PHI  PRIME 


DIMENSION  T(3l 
-  -  - .TAUI 


Bia§Rlf8R 


_  i( 19»90» 

DIMENSION  ISUMMPI 191 .ISUMMNI 19). I  ROW  I  50) 
equivalence  ((CCT.KOUNTY)  .(FHI.F)  , 


(ICCSi 


I  .ICTETAO.CTO) 


1  000 

\m 

1003 

1007 

1008 

1010 

1011 

1012 

1013 

18}^ 

181! 


..(COUNT)  I  _ 

l^rsfETAOVSTO) .(CTETAP.CfP) .ISTEfAP.STP) .(CTtTA.CT). ISTETA.ST) .(SPH 
21 .SP) , I CPHI .CP) ,( CPHIP. CPP) .ISPHIP.SPP) 

FORMAT  (1216) 

'1. 


(6F12,0) 

'  (6El2.e) 


NWYP«I6i 


i2HC115X.2HC2) 


3 

4 


format 

FORMAT 

format  (2E12.8) 

96IO) 

FORMAT! IH  4X.3HLGB5X.3HLEN5X.3HLEI5X,3HMCP6X.2HLS5X.2HLD6X.5HLOSTA) 

^^ORMAT  (6H1  NO«I6.5H  NP*I6.6H  NWY«I6»6M  NWZ«I6.7H 
17H  NWZP-I6.8H  NWYZM*I6,8H  NWY2P»I6.7H  NEGZ-I6) 

FORMAT! IH  7X . 2HV014X iZHPHlSX . IHH 15X , 2HA 1 1 5X  .2HA2 15X  .< 

FORMAT! lH^7xl2(Hio{4X.4HYLlMllX,3HYSClZX,4HZL IM 1 1 X . 3HZSC 15X . IHB lOX » 
14HJUMP) 

FORMAT!5E17,8) 

F0RMAT!7E16.8) 

^F0RMaI^”''(40H0  HISTOCiRAM  FOR  POSITIVE  VALUES  OF  Y  ) 


00002 


00003 

mn 

00006 

00007 

00008 

888®' 

888 

888 

888 

00019 

00020 

00021 

00022 

00023 

00D24 

888^^ 

888IS 


00031 

00032 

00033 

00034 

88811 

00037 


888S? 

00042 

00044 

00046 

00047 
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•-•orvo  nrv^ooorv^r>r%r»r> 


i02<; 

1023 


FORMArdH  .1916) 

FORMAT  (40H1  HI5TC0RAM  FOR  NEGATIVE  VALUES  OF  Y 
READ  INPUT  TAPE  5. 1000. KCT .JUMP  .NET. NT 

READ  INPUT  TAPE  5 . IOC  1 . YL I M . YSC.2L IM.Z SC .F LAGM .F LAGP 
READ  INPUT  TAPE  5 . 1 001 . YPSC .ZPSC . YQR.ZOR 
READ  INPUT  TAPE  5 . 1 00 1 . S TO .C TO.DSS .DC P .B 
READ  INPUT  TAPE  5 . 1 00 1 . PO .ZO .PH.H .PBT . ZBT 
READ  INPUT  TAPE  5 . 1 003 . ( ZH I  1) .PL  I  I  )  . I  *  1  .NE T  ) 

READ  INPUT  TAPE  5.  1001.  XI. X2 
DUMMYxINTRMFI XI  .X2) 

KOUNTxO 

REWIND  27 
RYPSC-l.Q/VPSC 
BZPSC»I.O/ZPSC 
TNUMxl 5./6. )**5 
RCRT  xi«/6. 


00049 

00050 

00051 

00052 

00053 

00054 

00055 

00055 

00057 

00058 

00060 

00061 


00052 

00063 

00054 

00065 

00068 

00069 

00070 

m 

*00073 


**«**************»*****»*START  OF  A  NEW  CASE ******************** *000 74 
01  READ  INPUT  TAPE  5 . 1002 . VO. A 1 . A2 .Cl .C2 .F 


NE  =  0 

00077 

NOxH 

00070 

N0*0 

00079 

Lr,B  =  o 

00080 

NxO 

00081 

LEI=0 

00082 

LEN  =  0 

00083 

LOSTAxO 

00084 

MCP  =  0 

00085 

LS  =  0 

00086 

LD  X  0 

NWZ  =  0 

00087 

NWY  =  0 

00088 

NCGZxO 

00089 

NWYZMxO 

00090 

NWYZP*0 

00091 

NWZP*0 

00092 

NWYPxQ 

00093 

*********************************************************** 

******** ********0009 5 

*#******•**#*»*#**»******■»»********•»***»••**#******** ********************00096 

10  N*N+1  00097 

IC*0  00098 

L*0  00099 

XN«0  00100 
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YN*0 
ZN  =  P 
PN  =  0 
VN  =  0 
00  =  11 
Vd  =  0 
VH*0 
XNB  =  0 
vnb*o 
7N9=ZO 
VNB*V0 
PNa=Po 

PARTICLES  BEEN  PROCESSED^ 


IF(N-NT 101 IStOl 15.0150 


00101 
00102 
00103 
OOlOA 
00105 
00106 
00107 
00108 
00109 
00110 
00111 
00112 
con  3 

»*0011A 

83HI 


>•«»**  It 

.0  continue  **  **  00120 

DO  152  IV  =  1,1=) 
isummp( iY)=c 

ISUMMNI IY)=o 
00  152  tz=5,AF 

I iUMMP (IV)  :WP (  I Y  ,  1 Z  )  *  I SUVMP I  I Y ( 

I  SUMMN(  tVlsNNI  IY,I2)  +  I  SU'AMN)  I  Y  ) 

152  CONTINUE 

DO  155  12=1,50 
IROw( IZ|»0 
DO  155  IV=1,19 

I  RO>v(  IZ  )  =  IROW(  U  )+VIP  (IV,  U  l+MNI  I  Y,  I  Z  » 

155  continue 

*RirE  OUTPUT  TAPE  6,1010,NC,NP,NwY,N»Z,N'a'VP,NwZP,N(,YZM,N'aVZP,NEGZ  00121 

write  output  tape  6.1008  00122 

WRITE  Output  tape  6, 1007  ,LGct,LtN,LE  I  .''CP  ,LS  .LD.LOSTA  00123 

WRITt^  OUTPUT  TAPE  6,1011  0012A 

WRITE  OUTPUT  tape  6 , 1 0  1 5 , VO ,PH ,H  ,  A 1 , A 2 , p 1 , C 2 

,  WRITE  OUTPUT  TAPE  6,1013  00126 

write  ouTPuT  tape  6, 1012 ,Z0,YLIM,Y5C,Zl Im.ZSC  ,3, JUMP 

WRITE  OUTPUT  TAPE  6,1021  00128 

WRITE  OUTPUT  TAPE  6,1C22,MP  00129 

WRITE  OUTPUT  TAPE  6,1023  00130 

WRITE  OUTPUT  TAPE  6,1022,MN  00131 


write  OUTPUT  TAPE  6,  lOK 

WRITE  OUTPUT  TAPE  6,  1022 

WRITE  OUTPUT  TAPE  6,  1022 

WRITE  OUTPUT  TAPE  6,  1022 

no  AOi^n  1  =  1,10 
DO  flu  .0  J=1 ,50 
MN ( I , J ) =0 
MP( I , J) =0 
CONTINUE 
END  file  27 


00121 

00122 

00123 

0012A 


00128 

00129 

00130 

00131 


I SUMMP 
1 SUMMN 
I  ROW 


00133 

00134 

00135 

00136 
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YM  =  0 
ZN  =  0 
PN  =  0 
VM  =  0 
noso 
Vd  =  0 
VH»0 
XN3=0 
VNB*0 
2NB=20 
VNRsVO 
PNBxPO 
r 

C***********»»***hAVE  NT  PARTICLES  BEEN  PROCESSED* 
r 

IF(N-NT)0115t0115.0150 


00101 

00102 

00103 

OOlOA 

00105 

00106 

00107 

00108 

00109 

00110 

00111 

00112 

com 

*0011A 

8SH^ 


iO  continue 

DO  152  IY=1,13 
isu*^vp(  iY)=c 
ISUMMNI IY)=^ 

00  152  12*5. A'' 

ISUMMP  (  I  Y  I  =^^P  (  I  Y  ,  1  2  )  *  I  SUMMP  (  I  Y  ) 

ISUMMNI  IY|=mniiy,I2)  +  I  5:jmmn(  I  y» 

152  continue 

DO  155  12=1 .50 
!RO«( 12  )«0 
DO  155  IY=1.19 

IrtO-v(  12  )  =  IROW(  12  I  +MP(  1  Y,  12  I't-MNI  I  Y.  12  I 
155  CONTINU'E 

aRITE  output  tape  6.1010.NC.NP»N»iiY»N»2.N'A’YP,Nw2P.M»Y2M,N'AYZP,NEG2 
aRITE  output  tape  6.1009 

aRITE  output  tape  6. lOOT.LGd.LEN.LE 1 .«CP .L'.LD.LOSTA 
WRIT*-  OUTPUT  TAPE  6.1011 
ARITE  OUTPUT  tape  6. 1015. VO. PH, h.AI ,A2,ri ,C2 
WRITE  OUTPUT  TAPE  6.1013 

aRITE  OUTPUT  TAPE  6 . 1 0 1 2 . 20 » YL I M . Y5C . 2 L I M , 2 SC . 3 . JUMP 
WRITE  OUTPUT  TAPE  6.1021 
aRITE  OUTPUT  TAPE  6.1C22.MP 
WRITE  OUTPUT  TAPE  6.1023 
WRITE  OUTPUT  TAPE  6,1022.MN 
aRITE  OUTPUT  TAPE  6.  1016 
WRITE  OUTPUT  TAPE  6.  1022.  I SUMMP 
WRITE  OUTPUT  TAPE  6.  1022.  I SUMMN 
write  OUTPUT  TAPE  6.  1022.  I  ROW 
no  snnn  1  =  1, IR 
DO  fli"'  .0  J=  1 , 50 
MN (  I  , J  )  =0 
MP( I ,J) =0 

JOCO  CONTINUE 

END  EIlE  27 


.**♦*******00118 

.**********00119 

00120 


00121 

00122 

00123 

0012A 


00120 

00129 

00130 

00131 


00133 

0013A 

00135 

00136 
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KOUNT^KOUNT+l 

IF  (K0UNT-<0UNTY)101.151.151 
0151  JUMPr2 

GO  TO  0101 


mn 


C  00142 

C  00143 

C  00144 

C  00145 

C  00146 

C  00147 

E  881tl 

C  00150 

C  00151 

C  00152 

C»*«*******************calCULATION  of  decay  TIMES******************00154 
r* •••*••**•••••••*••00155 

0115  RN1=URNDMF(  1.0)  00158 

RN2  =  URNDMF(  1.0)  00159 

RN3=URN0MF( 1.0)  00160 

1F(RN1-RN2)30.0115.40  00161 

30  IF(RN2-RN3)3l.0115.32  00162 

31  R11)=RN1  00163 

R( 2 ) =RN2 

R(3)=RN3  00165 

GO  TO  45  00166 

32  IF(RNl-RN3)33.0U5t34  00167 

33  R(1)=RN1  00168 

R(2)=RN3  00169 

RI3)=RN2  001.70 

GO  TO  45  00171 

34  Rll)=RN3  00172 

R(2)=RN1  00173 

R(3)=RN2  00174 

GO  to  45  00175 

40  IF (RN2-RN3)42.01 15.41  00176 

41  R(1)=RN3  00177 

R(2)=RN2  00178 

R(3)=RN1  00179 

GO  to  45  00180 

42  IF (RN1-RN3 )44.nl 15,43  Q0181 

gfJl®  88111 

R(3)=RN1  00184 

GO  to  45  00185 

88ii^ 

R(3)=RN3  00188 

45  DO  0117  1=1,3  00189 

Rl^sp  (  I  ) 

IFIRR-RCRT)  50.50,51 

50  TIME) I )=6.*RR 

GO  TO  0117  00193 

51  RR=1.-RR 
R2=RR*RR 
R5»R2*R2*RR 
TIME! I )=TNUM/R5 

0117  continue  00196 
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RTD=T IME ( 1 1 

00197 

TDTWO=  T IMc (2 1-T IME( 1 1 

00198 

T0ThR=TImE(3 )-TlME(2) 

00199 

0  1  ?c 

CONT 1 NUF 

00200 

*****»v»*rH005£  INITIAL  DIRECTION  OF  MOTION  AND 

TYPE  OF  PART 1CL£**********00201 

0113 

R=URNDMF (1,0) 

00202 

IF(R-F)P121 .0118,0131 

00203 

0121 

I  =0 

00204 

0  12? 

RrURND''^'  (1.0) 

00205 

CTETA=R 

00206 

STFTA=50RTF( 1.0-R**2) 

00207 

C  125 

R1=URNDMF(  l.(^) 

00208 

R2=URNDMF( 1.0) 

00209 

T  =  R  1**2  +  R2**2 

00210 

IF{T-1.01 0126.0126,0125 

00211 

0  126 

rphl=(Ri**2-R2**2)/T 

00212 

SPHI =2.C*R1*R2/T 

00213 

0127 

R=URND^'F  (1,0) 

00214 

IF(5-C.5)0129,'^l?').012e 

00215 

0128 

CiPH  I  i-SPMI 

00216 

0129 

FjO  =  N'^  +  l 

00217 

GO  TO  (0?'^1  ,66''  1  )  .JUMP 

00218 

0131 

013^ 


GO  TO  0'»ni 


1  =  1 

0:UPND“F  (1,0 
CTtTAP=R 

STC  TAPr'.QPT'^  (  1 ,0-R**2) 

,'.P  =  NP*  1 


C«***»i»CALCULAT  ION  OF  TIMES  TO  COLLISION 

0201  RfJ  =  URNDMF  (1.0) 

PN  =  PfKlB  +  Al*CTFTA*LOGF(RN) 
r  »**»**»*»**»-♦’»  does  PARTICLF  ESCAPE  WITHOUT 
I  F ( PN-PH)  0601  .OSOl .POA 
20A  IF(Pn-pbT)  2n2.il0,iir 
2C2  ASSIGN  203  TO  IH 

GO  TO  5600 

203  DZO=ZN-ZNB 

DSO=DZO/CTeTA 

0208  DXO=DSO»STf T A»CPHl 


00219 
00220 
00221 
00222 
00223 
0022A 
00??5 
00225 
00227 
00228 
00229 
00230 
00231 
00232 
00233 
»00234 

CF  NEUTRAL  PAR T 1  CL E S«*****»* *00235 
♦,»»*»»*»«*»****«»**«**»****»*00236 

00237 

COLLI S I  ON* ••*•**♦** ********** *002 39 

00240 

00241 

00242 

00243 

00244 

00245 

00246 
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DYOaOSO*STETA«SPMl 
yN  =  YN9+r)Y0 
XN  =  XNB4-0X0 
TsVNB**2-*'Cl*L0GF  (RN) 

IF  ( T  1021 At021A.n??0 
0214  NE*NF+3-L 

L0STA=1.0STA  +  1 

C***«**PaRTICLE  slowed  down  to  rest. end  history. si 
GO  TO  0110 

0220  VN*SORTF(T) 

VB=0.5*( VN+VNB) 

RTC«DSO/VB 

GO  TO  0401 

0250  WRITE  OUTPUT  TAPE  6.0251.K.L.I 

0251  FORMATOIIO) 

GALL  nu«P 


00247 

00248 

00249 

00251 

00253 

MIlARLY  at  OTHER  P0INTS**00254 

00255 
00256 
00257 
00258 
00259 
00260 
00261 
00262 


C«************»**********««****« »•*»»»»»•*•»*•«*  »  *»***«»*****************##5o253 

C*******i****«*»»*cALCULAT ION  OF  DECAYS  OF  NEUTRAL  PART  I CLES******»*********00264 
d  ***********************  *»«»*ii*»»*»i»*»«»*»»»*»***(n(****»***********#*******00265 
0401  CONTINUE  00266 

IF  (RTD-RTCI0440. 0440. 0403  00267 

d*»***********.«****».*****c0LLlS10N  BtFCRc  0cCAY***«»»*»*«*»*»»******<*****00268 
0403  . 


0415  LGB»LGB+1 

GO  TO  on 

C;»*«»**»ENo  MIS 
0420 


0440 


439 

441 


XNB  =  XN 

00269 

YN8*YN 

00270 

ZNB=2N 

00271 

PNP«ON 

00272 

VNB*VN 

O0273 

00274 

00275 

RTD»RTD-RTC 

00276 

CTETAPsCTETA*CTETAO-STETA*SPHI *STETA0 

00277 

15.0415,0420 

00278 

00279 

00260 

PROGRAM** *********** 002 8 I 

TORY,  similarly  at  all  other  points  of 

STETAP=SQRTF ( 1 . 0-C TE T AP»»2 ) 

00282 

1 

00283 

*»»*»»*««**»DFCAY  BEFOrtt  COLLI S I0N»**»*< 

NF=NF+1 

00286 

L^L+l 

00286 

i  =  m 

00287 

DSO=RTD*V° 

00288 

OZO=DSC*CTETA 

00289 

DXOsDSO*STETA*CPMl 

00290 

DYO=DSO*STET A»SPHl 

00291 

XNsXNB+DXO 

00292 

YNsYNB+OYO 

00293 

ZN*ZN3*OZO 

00294 

)  110.110.439 

00295 

1  TO  IB 

00296 

5500 

0P  =  DNj_PM9 

00298 

T=VNB**?+B»DP/CTETA 

00299 

62 


OAftO 

0470 

1 

0486 

0486 


0490 

0491 

0495 

0496 

0301 

f  •  **» 

0307 

3C8 

310 


0314 

0320 


IF  ( T 10460.046^.0470 
NF»NE*3-L 
LOSTAsLOSTA+I 
GO  TO  0110 

VN=5QRTF(T) 

V'3s0.5*  ( VN+VNB  1 
ASSIGN  1  TO  MCALL 
GO  TO 

continue 

IF  (L-3)OaP5. 04(16,0250 
LO=LO*1 

GO  T'l  0110 

XNO=XN 
YN3=YN 
2Ne  =  2:'l 
PNa=PN 
VNO=VN 

CTE  TAP  =  C  TET  A*CTcTAO-STETA*':PmI  * 
I F(CTP) 0415.0415.0490 

ST£TAP=SQfiTF ( 1 .0-CTE TAP*»2 ) 
RTC^PTC-RTO 

IF  (1-2)0495.0496.0250 
PTC=T0TW0 

GO  TO  0501 

RT0=T0Th9 

***»I*«2*  *#»**«.«  4  *«(.*»*«»••«»«  ««*•»<*#»** 
**4»»»CALCULATI0N  of  times  to  collision  of 

RNsijRNDvFI  1,0) 

DP=A2*CT£TAP*CT£TAO*LOGF(RN) 

PN»PN9*DP 

*4**»*****  does  PARTICLE  ESCAPE  WITHOUT  CC 
IF (PN-PH) 0001.0801 .0507 
IF(Pn-P9T ) 308.110, 1 10 
ASSIGN  310  TC  IH 
GO  TO  5600 

DZP=(ZN-2N9)/CTtTAO 

YN=YNB-DZP*STETAO 

XN=XN6 

T=VNB**2+C2*L0GF(RN) 

IF  ( T 10314.0314 .0320 
NE=NF+3-L 
L0STA  =  L0STA-.1 
GO  TO  0110 

VN=50RTF ( T ) 

VP=0,5*(VN+VN5) 

nSP=OZP/CTFTAP 

RTC=DSP/V9 

GO  TO  0501 


00300 

003of 

00302 

00303 

00304 

00305 

00306 

00307 

00308 

00309 

00310 

00311 

00312 

00313 

00314 

00315 

00316 

STETAO  00317 

QQ318 

00319 

00320 

mi 

00323 

003.24 

00325 

*»*»»»» »»**»»***4*«#********00326 
lOMSEO  PARICLtS«**»*****»*00327 
«»»»*»*«******«**»**********0(5328 

00329 

00331 

LLI  r-  1 0N**» •*«**«*****»» *****00332 

00333 

00334 

00335 

00336 

00337 

00338 

00339 

00341 

00342 

00343 

00344 

00345 

00346 

00347 

00348 

00349 

00350 

00351 

00352 

00353 

00354 
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00355 

00356 

mn 

00359 

00360 

00361 

*88161 

*00364 

»00365 

00366 

00367 

00366 

*00369 

00370 

00371 

00372 

00373 

00374 

00375 

00376 

00377 

00378 

00379 

00380 


►CALCULATION  OF  OECAVS  FOR  IONIZED  PARTICLES* 


0501  CONTINUE 

IF  (RT0-RTC)O540. 0540. 0503 
r**«***»***«»************»*cOLLI SION  BEFORE  DECAY* 
503  XNB=XN 

YNB=VN 
ZNB»ZN 
PNBaPN 
VNB«VN 
X«K+1 
I«l-I 

RTOsRTD-RTC 
IF  (  I  10250.0510.0301 

RNAaURNOMFI 1 .0) 

RNB=URNDMF(1.0) 


0510 


0»RNA**2-*PNB**2  00381 

IF(Q-1. 010515. 0515. 0510  00382 

0515  CPHIPs(RNA**2-RN8**2 » /Q  00383 

SPHIP52.0*R\A*RNB/Q  00384 

0516  R»URN0MF( 1.0 1  00385 

IF(R-0. 510517. 0516. 0518  5038§ 

0517  SPHIP.-SPHIP  00387 

0518  CTETA*CTETAP*CTETA0*STETAP»SPHIP*STETA0  00388 

IF(CT)0520. 0520. 0525  92252 

0520  LGB»LGB+1  2222? 

GO  TO  one  ,  9922i 

0525  STFTA«SORTF( 1.0-CTFTA**2)  00392 

0528  R»URN0MF( 1.0)  2922? 

CPHI«STETAP*CPHIP/STETA  00394 

SPHI.SORTFl 1.0-CPMI**2)  2922? 

IF(R-0. 5)0530. 0528. 0531  29225 

0530  SPHI.-SPril  90221 

0531  GO  TO  I  0201 .6601 )  .JUMP  00398 

X* ********************* **»»0EC AY  bEFORc  COLL  I S ION* ********************* ****00399 
0540  NE*N|n  0040J 

I»m  00402 

0SP*RTD*VB  22^2? 

D2P=OSP*CTETAP*CTETAO  00404 

ZN-ZNB+DZP  22t2? 

IF(ZN-ZBT)110.11C.549  22?25 

549  IF(ZN-H)0541 .0550.0550  29^2Z 

550  CONTINUE  00408 
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rv^  r>  o  o  or»  o  o  rsr»  o  ri 


NE*NE-1 

00409 

U*L-1 

00410 

1  =  1-1 

00411 

DP=PH-PN8 

00412 

T=VNB**2+B*DP/(CTETAP»CT£TA0» 

00413 

IF(T)0314. 0314, 0805 

00414 

0541 

XN=XNB 

00415 

YN»YNB-DSP*CTFTAP*STETA0 

00416 

42 

ASSIGN  543  TO  IP 

00417 

43 

GO  TO  5500 

DPsPN-PNB 

mil 

TsVNB**2+(B*DP) / (CTETAP*CT£TAO) 

00420 

IF 

11)0560,0560,0570 

00421 

0560 

NE=NE+3-L 

00422 

losta«losta*i 

00423 

GO 

TO  0110 

00424 

057C 

VN»SORTF(T ) 

00425 

y*l’0.5*(  VN*VNB  ) 

ASSIGN  2  TO  MCALL 

TS 

GO  TO  7''00 

00428 

2 

continue 

00429 

0586 

IF(L-3)0585. 0586, 0250 

L0*L0+1 

00430 

00431 

GO 

o 

c 

o 

00432 

0585 

XNB*XN 

00433 

YN6»YN 

0043'4 

2N0»2N 

PNB=PN 

00435 

00436 

RT?  =  R1*C-«T0 

0Q437 

00438 

IF 

00439 

00440 

0595 

GO 

TO  0501 

00441 

0596 

rtd=tdthr 

00442 

GO 

TO  0501 

00443 

00444 

00445 

00446 

00447 

00448 

00449 

00450 

00451 

«*«»#*»#*  «»»********»***»*»»  #QQ452 

n.****!*  j$I*2*in»***00454 


mil 

6601 

CONTINUE 

00457 

DP=PH-PN0 

00458 

DZ=M-ZN3 

00459 

DS=DZ/CTETA 

00460 

dx=os*steta*cphi 

00461 
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DY=DS*STETA*SPHl 

Q0^62 

XN=XNB+DX 

00463 

YNsYNB+DY 

00464 

T*VNB**2+B»DP/CTETA 

00465 

lFrT)6250. 6250. 6470 

00466 

III? 

IF(2N)0260.6456.0460 

2N=H 

mil 

GO  TO  0460 

00469 

6470 

VH=SQRTF(T) 

VB=0.5*(VH+VNBI 

RTE=DS/VB 

00470 

00471 

00472 

00473 

6327 

I F(RTD-RTE >6440.6440.6450 

6440 

RTC=RTE 

00474 

GO  TO  0440 

00475 

6450 

RTCH=RTF 

00476 

GO  TO  0639  03477 


C************  ESCAPf  KOUTINE  FOR  NEUTRALF  ***»*»»•*•»*«•**«*»««*•**•***«  00479 

0*«*******»****«*«****Ov*»**»»«»««»*»»**‘'***><>»**«**«***>«*»*««********»«*C0480 


0601 

CONTINUE 

00481 

0P=Ph-PNB 

00482 

D2=H-ZNB 

00483 

DS=DZ/CTETA 

00484 

DX=DS*STETA*CPHI 

00485 

DY=05*STETA*SPHl 

00486 

XNsXNB+i^X 

03487 

YN  =  yn3-.0y 

00488 

TsvN9**2+9*OD/CTETA 

00489 

IF(T)0615. 0615. 0620 

00490 

0615 

IF(ZNI025?. 0616. 0214 

0049-1 

0616 

2N*H 

00492 

GO  TO  0214 

00493 

0620 

VH=SQRTF ( T ) 

00494 

Ves0.5*( VH.VNB) 

00495 

RTCH.DS/VB 

00496 

GO  TO  0630 

00497 

0630 

IF (RTD-RTCHl 0635 .0635.0639 

00498 

0635 

RTC*RTCM+9000.0/ (CTETA*VH) 

00499 

GO  TO  044'7 

00500 

39 

LEN»LEN+1 

00501 

0640 

RTS*OSS/(CTETA*VH) 

00502 

RTD=RTO-RTCH 

00503 

IF  (RTD-9TS)0650 .0650.0680 

00504 

0650 

NE»NE-M 

00505 

L«L-M 

00506 

2NB*H 

00507 

DS=RTD*VH 

0Z*DS*CTETA 

88181 

DX=DS»STETA*CPHI 

00510 

DY.05*5TFTA*5OHI 

00511 

2N=ZM6+DZ 

00512 

XN*XNB*DX 

00513 

YN*YNB+OY 

00514 

ASSIGN  3  TO  MCALL 

00515 

GO  TO  7000 

00516 

66 


3  CONTINUE 

IF(L-3)0665i0666.0250 
0666  Ln=Ln+i 

GO  TO  0110 

0665  XNB*XN 

YNB=YN 
?NEsZN 

CTETAPsCTETA*CTETA0-3TETA*«=HI*6TETAC 

IF(CTP|OA15«0A15t06/0 

0670  STETAP=5QPTF( 1.0-CTETAP*»2) 

:7TCp=dcf/(CTEtap*.vh) 

0675 

GO  TO  OOCO 

0676  PT0=TDTHR 

GC  TO  0900 

0680  L5=L5+1 

A6C0  CONTINUE 

GO  TO  0110 

C**«**»»it****j*<*»****«e  It  **»»»«»*»»&»  £*»«»»*»»•»♦•*♦•♦■• 

FSCAPE  ■’CUTINe  F09  lOMZf.C-'  PARTICLE 

^«tttttt«itirititttiti>itititititititit»tf->#n'ttt!tit»*v:ti^itiiit-,rit<S''ttit«if>f>'!;if«{.'!r 

0801  CONTINUE 

L£I=LEI+1 
0P:PH-?NB 

T=VNe*'>2+o^DP  /  (CTtTAPitCTETAO) 


00517 

00518 

00519 

00520 

00521 

00522 

00523 

00524 

00525 

00526 

00527 

§8111 

00530 

00531 

00532 

00533 

00534 

00535 

:  *»««««<*«*••»**•  00537 

i»»..*****.*»******»#*oo538 

00539 

00540 

00641 

00542 


0  8  04 
1  “06 

0£C5 


0820 

C825 


826 

827 

0830 

1836 

0  940 


{F(  T  i; 
IFIZN 


,  ,  10804. 0804f 0605 

t F(ZN)025^  .  1806.03  14 
ZN  =  H 

10  TO  P914 

VHrSORTF ( T I 
VP=0.5»( VH+VNB I 
CZ=H-ZNB 

DS  =  DZ/(CTE  TAP*CTETAO) 
RTH=D5/V5 


GO  TO  0820 

IF  (RT0-RTh)C325. 0825. 0870 


NE  =  N‘^  +  1 
L  =  L+1 

DSP=PTD*VP 

0ZP=0SD*CTETAP*C7ETA0 

XN=XNB 

YMrYNB-DSP*CTPTAO»6TFTAO 

ZN=ZNH+DZP 

IF(ZN-ZRT)  110.110.826 
ASSIGN  827  TO  IP 
GO  TO  5500 

DP=PN-PNB 

T=VNB*»2+(P«DP )/ (CTtTAP*CTETAOI 
IFIT)0830. 0830. 0840 
IF(ZN)025C. 1836. 0560 
ZN  =  H 


GO  TO  0560 


VN=SORTF ( T  ) 


§§! 


3543 

5544 

00545 

00546 

00547 

00540 

00549 

00550 

00551 

00552 

00553 

00554 

00555 

00556 

00557 

00558 

00559 

00560 

00561 

00562 

00563 

00564 

00565 

00566 

00567 

00568 

00569 

00570 
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Ve=0.5*( VN+VNB) 

ASSIGN  A  TO  mcALL 
GO  TO  7000 
CONTINUE 

IF(L-3)0850.08S1  .0250 
LD»L0+1 

GO  TO  0110 

XNBsXN 
YNB=YN 
ZNB=2N 
PNB=PN 
VNB=VN 
RTH=RTH-RTD 
IF  (1-2)0860.0862.0250 
RTO=TOTWO 

GO  TO  0820  06*36 

RT0=TDTHR  00387 

GO  TO  0820  00588 

RTC*RTD-RTH  00589 

RTCP=DCP/(CTETAP*VH)  00590 

2NB=H  00591 

DSP=RTH*V6  00592 

ynB»YNB-DSP*CTETAP*STETAO  00593 

GO  TO  0900  00594 

IF  (RTD-RTCP)09C5. 0905. 0950  00595 

NE'NE'.l  00596 

L*L+1  00597 

0SP=RTD*VH  00598 

D2P»OSP*CTETAP*CTETAO  00599 

XN*XNB  00600 

YNsYN8-DSP*CT£TAP*STfcTAO  00601 

2N«2NB+D2P  00602 

ASSIGN  5  TO  MCALL  00603 

GO  TO  7000  00604 

CONTINUE  00605 

IF(L-3)0920. 0919.0250  00606 

LD«LO*l  00607 

GO  TO  0110  00608 

RTCP*RTCP-RTD  00609 

XN0=XN  00610 

r.ur.  88111 

IF  (L-2)0925. 0926. 0250  00613 

RTO=TDTWO  00614 

GO  TO  0900  00615 

RTD=TOTHR  00616 

GO  TO  0900  00617 

continue  00618 

VCP*MCP+1  00619 

CONTINUE  00620 

GO  TO  0110  00621 

00622 

****************************************************************** ****00624 


00571 

00572 

00573 

00574 

00575 

00576 

00577 

00578 

00579 

88ti? 

00582 

00583 

00584 


C 

c**** 


IF 

(L- 

■2)0925.0926.0250 

0925 

RTO=TDTWO 

GO 

TO 

0900 

0926 

RTD=TOTHR 

GO 

TO 

0900 

950 

continue 

VCP^MCP-H 

4900 

CONTINUE 

GO 

TO 

0110 
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c  ************  INTERPPOLAT ION  SUBROUTINE  FOR  PRESSSURE*****«**** 

5500  IF  (ZN-ZH(in  552O.551Ot5501 
55C1  IF (ZN-ZH( NET) >5502.5515 .5523 


5502 

5503 

5  5  04 
5510 

5515 


5520 

5521 

5525 

5526 


DO  5504  1*1. NET 


IF  iZN-ZHIIl)  5503.5504.5504 
PN=( (ZN-ZH( i-1 ) )*PlTI I-TZN-ZM 


( I  I )*PL1 l-l )  ) / ( ZH(  I  )-ZH( I-l )  ) 

PN=10.**PN 
GO  TO  5599 
CONTINUE 

PN*10.**PN 
GO  TO  5599 
PN=PL(NET) 

PN*10.**PN 
GO  TO  5599 

WRITE  OUTPUT  TAPE  6.5521.ZN 

FORMAT! IH  3X.9HALTTTU0E=E16.8.21H  BELOW  RANGE  OF  TABLE) 

GO  TC  5510 

WRITE  OUTPUT  TAPE  6.5526.ZN 

-  .  .  . 


FORMAT! IH  3X.9HALTITUOE 
- 551T 


628 
00629 

8811? 

00632 

00633 

00634 

00637 

00638 

00639 

00640 

00641 

00642 

00643 

00644 

00645 

00646 

00647 


_iE16.8.21H  ABOVE  RANGE  OF  TABLE) 

GO  TO  5515 

5599  GO  TO  IP.  1441.543.827) 

r****»***«**************#*********4******«*»*»****«******** *•**************£ 

C  »*»»».»>•»••  INTERPPOLAT ION  SUBROUTINE  FOR  ALT  I TU0ES****«*********00650 


5600 

5601 

5602 

5603 


5604 

5605 
5610 
5615 

5620 

5621 

5625 

5626 

5699 

C 

c 

c 

c 

c 


X*L0GF!PN)*. 43429448 

IF!X-PL!1))  5602.5610.5620 
IF  !X-rLINET))  5625.5615.5603 
OC  5605  1*1. NET 

IF!X-PL!I))  5605.5605.5604 

ZN*!  !X-PLI I-l ) )*ZH! I )-!X-PL! I ) )*ZHl I-l ) I / IPLI I  )-PL! I-l) ) 
GO  TO  5699 
CONTINUE 
GO  TO  5699 


00648 

*00649 

*00650 

‘88^1^ 


00653 
00654 
00655 
00656 
00657 
00658 
00659 
00660 
00661 
00662 
00663 
00664 
00665 
00666 
00667 
00668 
00669 
00670 
00671 
00672 
00673 
00674 
00675 
00676 
00677 

r****»***»******»»****»*******»»***»***»*»»»***»*****»*************«*******00678 
r  *»**»»»»******PREPaRAT I  ON  OF  OUTPUT  maT R I X ****** ****♦♦**♦**♦  00679 
C* *«**************»***«****»******»**«*•»•***<«*•*******************•****** 00680 


ZN*ZH! 1 ) 

GO  TO  5699 
ZN*2H!NET) 

GO  TC  5699 

WRITE  OUTPUT  TAPE  6. 5621. X 

FORMAT! IH  3X.6HL0G  P=E16.8.21H  ABOVE  RANGE  OF  TABLE) 
GO  TO  5610 

WRITE  OUTPUT  TAPE  6. 5626. X 

FORMAT!lh  3X.6HL0G  P=tl6.8»21H  BELOW  RANGE  OF  TABLE) 
GO  TO  5615 
GO  TO  IH.!203.310) 
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7  0?0 
7C09 


7010 

7008 

7011 

7019 


7027 


7026 

7C21 

7C22 

7120 


CONTINUE 

IF  (FLAGP)  7009. 7020. 7010 
YP»( YN*CT0I+( t2N-ZO>«STO) 
2P=1 (ZN-20)*CTO)-( YN*STO) 
GO  TO  7008 
YP  =  YN 


7050.7050.7011 

7060.7060.7012 


IF?RY^§C  -ABSF I YP) I 
IFIRZPSC-ABSF(ZP) I 
CONTI NUF 

IF (FLAGMl 7019.7999.7027 
YP»YN 
ZP*ZN 

GO  TO  7021 

continue 

ZP*(  tZN-ZOI*<.TO)-l  YN*STO) 

YP»( YN*CT0)+( (ZN-201*ST0I 
IFIZPI  7026.7026.7021 
NEGZ=NEG2>1 
GO  TO  7<J99 
IZxIZP/ZSCl+l.O 

IFr49-I2)  7030.7022.7022 

IF ( 18-XABSF( I Y 1)7035.7125.7125 


00681 

mu 

00684 

00685 

00686 


00689 

00691 

00692 

00693 

00694 

00695 

00696 

00697 

QQ698 

00699 

00700 


7125 

7130 

7132 

7C3C 

7031 

7035 

7036 

7050 

7051 

7060 

7999 


» 


MN( lY.IZ I »MN( lY.IZ )+l 
GO  TO  7999 
1  Y*YP/YSC  +  l.g 

IF( 18-XA3SF(Ty 1 17035.7132.7132 
MP( Iy.IZ 1 »MP( I Y. IZ 1+1 
GO  TO  7999 
NWZ3NW2>1 

IFI 19-XABSF( lY 1 17031.7999.7999 
NWYZM.NWYZM+l 
GO  TO  7999 
NWY=NWV+1 

IF(50-IZ)  7036.7999.7999 
N»(YZM*NWYZM't-l 
GO  TO  7999 
NMYPzNWYP^I 

IFIRZPSC-ABSFIZP) 1  7051.7999.7999 
NWYZP-NWYZP.l 
GO  TO  7999  • 

NWZP*NWZP+1 
GO  TO  7999 
CONTINUE 

GO  TO  MCALL. ( 1.2.3.4.51 

END 

DATA 


00712 

00714 

mi 

00718 

00719 

00722 

00723 

00724 

00725 

00727 

00728 
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LMSC-B007043 


list  of  ^mbols 


KCT 

JUMP 

NET 

NT 

Y  LIM 
YSC 
Z  LIM 
ZSC 

FLAGM 

FLAGP 

YPSC 

ZPSC 

YOR 

ZOR 

STO 

CTO 

DSS 

DCP 

B 


PH 

H 

PBT 

ZBT 


Total  card  count 

Indicator  cross  section 

Number  of  entries  in  table  of  atmosphere 

Total  number  of  particles  to  be  processed 

Highest  value  of  Y*  in  output  matrix 

Size  of  grid  on  Y'  in  output  matrix 

Highest  value  of  Z'  in  output  matrix 

Size  of  gprid  on  Z'  in  output  matrix 

Indicator  for  matrix 

Indicator  for  plot  routine 

Y  scale  factor  for  plot  routine 
Z  scale  factor  for  plot  routine 

Y  position  of  origin  on  plot 
Z  position  of  origin  on  plot 
Sin 

Cos 

o 

Distance  beyond  which  particle  is  ignored 
Distance  to  conjugate  point 
See  Eq.  (4.8) 

Pressure  at  starting  point 
Starting  altitude 
Pressure  at  cut-off  altitude 
Cut-off  altitude 

Pressure  at  lower  cut-off  altitude 
Lower  cut-off  altitude 

Initialization  valves  for  the  random  number  generator  (IBM  only) 


% 
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Section  6 


ATMOSPHEWC  PARAMETERS  AS  A  FUNCTION  OF  ALTITUDE 
NEAR  SUNSPOT  MAXIMUM 


altitude 

(km) 

pressure 

(dyne/cm2) 

100 

1.74 

X 

10"^ 

120 

3.4 

X 

10'^ 

140 

1.04 

X 

10-2 

160 

5.1 

X 

10-2 

180 

3.1 

X 

10-2 

200 

1.95 

X 

CO 

1 

o 

220 

1.20 

X 

10-2 

240 

8.5 

X 

1 

o 

sH 

260 

6.4 

X 

1 

O 

280 

4.7 

X 

1 

o 

300 

3.6 

X 

10-“* 

320 

2.7 

X 

1 

O 

340 

2.04 

X 

1 

O 

s-« 

360 

1.54 

X 

O 

1 

380 

1.23 

X 

1 

o 

400 

9.8 

X 

10-2 

450 

5.2 

X 

10-2 

500 

2.9 

X 

10-2 

600 

1.00 

X 

10-2 

700 

3.5 

X 

10-2 

800 

1.32 

X 

10-2 

900 

4.9 

X 

IQ-'^ 

1000 

1.90 

X 

o 

1 

1200 

3.2 

X 

O 

1 

<x> 

1400 

6.7 

X 

10-2 

1600 

2.1 

X 

10-2 

1800 

1.14 

X 

1 

o 

2000 

9.5 

X 

10-12 

2500 

7.2 

X 

10-12 
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ATMOSPHERIC  PARAMETERS  AS  A  FUNCTION  OF  ALTITUDE 
NEAR  SUNSPOT  BflNIMUM 


altitude 

(km) 

pressure 
(dyne/cm* ) 

100 

1.74 

10-^ 

120 

2.1 

1 

O 

140 

4.6 

10-3 

160 

1.86 

10-3 

180 

9.1 

10-* 

200 

5.0 

10-* 

220 

2.8 

10-* 

240 

1.77 

10-* 

260 

1.14 

10-* 

280 

7.6 

10-3 

300 

5.1 

10-3 

320 

3.5 

10-3 

340 

2.34 

10'3 

360 

1.66 

10-3 

380 

1.14 

10-3 

400 

8.3 

10-3 

450 

3.6 

10-3 

500 

1.66 

10-3 

600 

3.4 

10-”^ 

700 

7.9 

10-3 

800 

2.4 

10-3 

900 

1.26 

10-3 

1000 

9.8 

10-3 

1200 

7.2 

10-3 

1400 

6.2 

10-3 

1600 

5.1 

10-3 

1800 

4.3 

10-3 

2000 

3.8 

10-3 

2500 

2.8 

10-3 
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